AIAA Journal 

Vol. 48, No. 7, July 2010 


Micromechanics Modeling of Composites Subjected to 
Multiaxial Progressive Damage in the Constituents 


Brett A. Bednarcyki 

NASA John H. Glenn Research Center at Lewis Field, Cleveland, Ohio 44135 

Jacob Aboudi- 

Tel-Aviv University, 69978 Ramat-Aviv, Israel 
and 

Steven M. Arnoldi 

NASA John H. Glenn Research Center at Lewis Field, Cleveland, Ohio 44135 

DPI: 10.2514/1.45671 

The high-fidelity generalized method of cells composite micromechanics model is extended to include constituent- 
scale progressive damage via a proposed damage model. The damage model assumes that all material nonlinearity is 
due to damage in the form of reduced stiffness, and it uses six scalar damage variables (three for tension and three for 
compression) to track the damage. Damage strains are introduced that account for interaction among the strain 
components and that also allow the development of the damage evolution equations based on the constituent material 
uniaxial stress-strain response. Local final-failure criteria are also proposed based on mode-specific strain energy 
release rates and total dissipated strain energy. The coupled micromechanics-damage model described herein is 
applied to a unidirectional E-glass/epoxy composite and a proprietary polymer matrix composite. Results illustrate 
the capability of the coupled model to capture the vastly different character of the monolithic (neat) resin matrix and 
the composite in response to far-field tension, compression, and shear loading. 


Nomenclature 


normalized instantaneous slope 


A, B 

= 

material parameters 

blb9 


normal tensile and compressive weighting 
factors 

D, H,L 

= 

dimensions of repeating unit cell, m 

Dj,Df 

= 

tensile and compressive damage variables 

da 

= 

increment in crack length, m 

da‘ M 


mode-specific increment in effective crack 
length, m 

d t 

= 

secondary damage variable 

dW‘ M 

= 

mode-specific increment in released strain 
energy density, Pa 

da, hfj, ly 

= 

dimensions of each subcell, m 

Ei.E 2 .E 3 

= 

Young’s moduli. Pa 

E? 

= 

current (damaged) Young’s modulus, Pa 

pO p 0 pO 

= 

initial Young’s moduli, Pa 

G 

= 

strain energy release rate, J /m 2 

G c m 


mode-specific critical strain energy release 
rate, J /m 2 

G' m 

= 

mode-specific strain energy release rate, J /m 2 

G 23 , G 13 , G 12 

= 

Shear moduli, Pa 

G c 

= 

Benzeggagh-Kenane failure criterion mixed- 
mode critical strain energy release rate, J/m 2 

y —0 /- 0 ^-0 

u 23’ °"13’ '~ r l2 

= 

initial shear moduli, Pa 


kj = instantaneous slope of material uniaxial stress- 

strain response, Pa 
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of material 




uniaxial stress-strain response 

li 

= 

material length, m 

Na, Np, Ny 

- 

number of subcells within repeating unit cell in 
each Cartesian coordinate direction 

Qe, S e 

= 

shear damage initiation strains 

Sit 

= 

compliance matrix terms, Pa~' 

S°j 

= 

initial compliance matrix terms, Pa -1 

t 

= 

depth of cracked body, m 

W 

= 

strain energy density released, Pa 

W. 

= 

total energy released, J 

W? 

= 

critical compressive strain energy, J 

xL xL yL 

= 

tensile and compressive nonnal damage 

YLZlZC 


initiation strains 

X T ,X C 

= 

axial tensile and compressive strengths, Pa 

Xl,X 2 ,X 3 

= 

global Cartesian coordinates, m 

yl.y 2 .y 3 

= 

local Cartesian coordinates, m 


= 

subcell local Cartesian coordinates, m 

a 

= 

power-law failure criterion exponent 

(afiy) 

= 

subcell indices 

Yu 

= 

engineering shear strain components 

AT 

= 

temperature deviation from reference, °C 

s ‘i 

= 

normal strain components 

C D C D C D 
b \ ’ fc 2 » fc 3 

= 

damage strains 

7 

= 

Benzeggagh-Kenane failure criterion exponent 

v 12 > v 21 > h’ 13 , 
v 31> v 23, v 32 

— 

Poisson’s ratios 

0 0 0 
y 12’ l 21> V 13 , 

v 31’ v 23’ V 23 

= 

initial Poisson’s ratios 

G ip T U 

= 

normal and shear stress components. Pa 

of 

= 

damage stress. Pa 


I. Introduction 

T HE importance of capturing the progressive nature of failure in 
polymer matrix composites was highlighted by the World-Wide 
Failure Exercise (WWFE) [JJ. While many ply-level failure criteria 
performed reasonably well when damage initiation was followed 
closely by final failure, in situations where significant nonlinearity 
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occurred before final failure, most standard criteria were insufficient. 
Methods that capture this nonlinearity have shown much better 
agreement with experimental data from the WWFE ( cf . , B ogetti et al . 
[2] and Nelson et al. [3]). 

In the present paper, the high-fidelity generalized method of cells 
(HFGMC) micromechanics model [4,5] is extended via an elastic 
progressive damage model to predict progressive failure of polymer 
matrix composites. This work was motivated by previous efforts [ 6 ] 
that indicated that HFGMC’s damage modeling capabilities were 
insufficient. The damage model employed previously within 
HFGMC relied on typical failure criteria (e.g., maximum stress, 
maximum strain, Tsai-Hill, and Tsai-Wu) to predict local failure 
within the subvolumes (subcells) within the HFGMC composite, 
repeating unit cell (RUC) based on the local (constituent level) 
stress-strain state. Then, upon failure, all of the subvolume’s stiffness 
components were instantaneously reduced to a very low value. The 
damage progression was thus modeled as a series of these instanta- 
neous subvolume failures before final failure of the composite or 
laminate. This previous work [ 6 ] revealed two major shortcomings of 
the approach: 

1) The damage evolution (essentially a step function) was not 
sufficient to capture the progressive nature of the failure, particularly 
in shear dominated cases. 

2) The single damage parameter that affected all stiffness 
components equally (i.e., lack of multiaxiality) was insufficient to 
capture the vastly different behavior observed in tension, compres- 
sion, and shear. 

To overcome these deficiencies, a damage model is proposed 
herein that accounts for the multiaxiality and progressive nature of 
damage within the constituent materials as simply as possible. It is 
thus assumed that all nonlinearity exhibited by the composite is due 
to damage, resulting in stiffness reduction, within the constituent 
materials. To account for differences in the response in tension, 
compression, and shear, three scalar damage variables and three 
damage strains are tracked in tension and three in compression. The 
effects of the damage variables on the material shear stiffness 
components can be controlled by additional material parameters. The 
damage strains are related to the orientation of the damage on the 
constituent scale and, assuming that these damage strains correctly 
capture the interaction effects among the strain components, it is 
possible to determine the evolution equations for the damage 
variables based on uniaxial stress-strain curves for the constituent 
materials. 

One way for a subvolume (subcell) within the composite to reach 
final failure in a given direction is for the damage to evolve until the 
subvolume stiffness in a given direction is close to zero. However, it 
is desirable to enable subcells to fail suddenly after a certain amount 
of damage has accumulated. As such, final tensile and shear failure 
criteria are introduced based on the mode-specific strain energy 
release rates, which can be calculated due to the directionality of the 
damage in the proposed model. For compression, a total dissipated 
strain energy criterion is considered. 

The proposed damage model uses the continuum damage 
approach first introduced by Kachanov [7,8], in which continuous 
damage variables are introduced to quantify the stiffness reduction of 
a material. The theoretical background of this approach, and 
specifically, the elastic damage approach, is discussed in detail by 
Lemaitre and Chaboche [9], Krajcinovic [J_0], Voyiadjis and Kattan 
[11], and Talreja [_12,_1_3]. One particular anisotropic elastic 
continuum damage model that has received a good deal of attention 
in recent years and is similar in some ways to that proposed herein is 
that of Matzenmiller et al. [141. This model has been incorporated 
within LS-DYNA [J_5,_16] as a progressive damage constitutive 
model for shell elements. The shell element progressive damage 
model available in ABAQUS ] 1 7 1 also uses an approach similar to 
Matzenmiller et al. [14], although the shear damage variable and 
damage evolution equations are different. 

Like the model of Matzenmiller et al. [14], the proposed damage 
model considers anisotropic elastic damage using a number of scalar 
damage parameters. Both models also track tensile and compressive 
damage separately. The model of Matzenmiller et al. [L4] was 


intended to simulate the response of an anisotropic composite ply, 
and thus considers plane stress conditions. The model proposed 
herein is intended for application to the constituents within a 
composite material, and thus considers a three-dimensional state of 
stress and strain. The damage initiation criteria employed by 
Matzenmiller et al. [ 14 ] are based on the Hashin [ 18 ] stress-based 
failure criteria, while the proposed model uses three-dimensional 
strain-based Hashin-like damage initiation criteria. Matzenmiller 
et al. [ 14 ] also include an independent in-plane shear damage 
parameter, while the model described herein relates the shear damage 
parameters to the normal damage parameters. Finally, the damage 
evolution equations and final-failure criteria of the two models are 
distinct. It is also noted that, like the model of Matzenmiller et al. 
[14], the proposed damage model could be readily implemented 
within finite-element models, including ABAQUS [ 17 ] and LS- 
DYNA [J_5], through their user material subroutines (although not 
done at this time). The present modular implementation of the model 
has been incorporated within NASA’s Micromechanics Analysis 
Code with Generalized Method of Cells (MAC/GMC) [19]. 

After the micromechanics theory and the proposed damage model 
are described, the coupled micromechanics-damage approach is 
applied to simulate a unidirectional F-glass/epoxy composite, for 
which the constituents are characterized based on ply-level data from 
the literature [1], The model’s ability to capture vastly different 
damage behavior in tension, compression, and shear is highlighted as 
monolithic and composite nonlinear stress-strain curves, and 
damage initiation and final-failure envelopes are presented. 

II. Micromechanics Model 

The HFGMC micromechanics model has been adopted to model 
the composite-level behavior in this paper. The HFGMC model is 
based on a homogenization technique for composites with periodic 
microstructure, as shown in Fig. la, in terms of the global coordinates 
(xj , jc 2 , X 3 ). The parallelepiped RUC (Fig. fb) defined with respect to 
the local coordinates (y, , y 2 , y 3 ) of such a composite is divided into 
N a ,Np, and N y subcells in the y 1; y 2 , and y 3 directions, respectively. 
Each subcell is labeled by the indices (cq0y) with a = 1, ... ,N a , 
yS = 1 , ... ,Np, and y= 1. ... ,N y and may contain a distinct 
homogeneous material. The dimensions of the subcell are denoted by 
d a ,hp, and / y , respectively. A brief outline of the HFGMC method is 
provided next. The reader is referred to Aboudi [4] for the full 
theoretical development in the case of linear electromagnetother- 
moelastic materials and Aboudi et al. [20,2j_] for the two- 
dimensional case of continuous fibers including inelasticity of the 
phases. 

1) A quadratic displacement field is assumed for each subcell in 
terms of the local subcell coordinates (y , y 2 ] , and y 3 y) in Fig. 1 ) and 
21 unknown coefficients. 

2) Equilibrium conditions, periodic boundary conditions, and 
interfacial continuity conditions of displacements and tractions 
between subcells are imposed in the average (integral) sense. 

3) The result is a system of linear algebraic equations, which is 
solved for the 2 1 unknown coefficients per subcell in terms of the 
constituent material properties, microstructural dimensions, and 
thermal and inelastic terms. 

4) Knowledge of these coefficients enables establishment of a 
localization relation that relates the average strain in each subcell to 
the externally applied average strain [5], along with thermal and 
inelastic terms. 

4) The local thermoinelastic constitutive equations are used in 
conjunction with the average stress relations to arrive at the effective 
thermoinelastic constitutive equations for the heterogeneous 
material. 

As stated, the original version of HFGMC requires solution of a 
system of 21 N a NpN y algebraic equations. However, just like the 
reformulation of generalized method of cells micromechanics model 
that has been presented by Pindera and Bednarcyk [22], it is possible 
in the case of perfect bonding between the phases to now use the 
continuity of displacements across the interfaces between the phases 
and significantly reduce the number of equations that have to be 
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c) 

Fig. 1 HFGMC model: a) a multiphase composite with triply periodic 
microstructures, b) the RUC, defined in the (y, coordinate 

system, is discretized into N a x Np x N r subcells, and c) the monolithic 
subcell is defined in the local coordinate system (y\ a> ,yf ’ ,y > i r> )- 


solved. This reformulation of HFGMC has been presented by Bansal 
and Pindera [ 23 ] in the special case of continuous fibers (i.e., the two- 
dimensional, doubly periodic case) and elastic phases, with the more 
general case presented by Arnold et al. [ 24 ] The unknowns in this 
reformulation are given by the average displacements at the surfaces 
of the subcells, which can be related to the original HFGMC 
unknown coefficients. However, by formulating the HFGMC theory 
in terms of the subcell surface average displacements, the size of the 
system of algebraic equations that must be solved becomes 


9N a N p N y + 3(N a Np + N a N y + N p N y ) 


rather than 2lN a NgN y , with an obvious significant savings. 
The ramifications of these savings are discussed in detail by Arnold 
et al. [24]. 

More recently, Haj-Ali and Aboudi [ 25 ] formulated the HFGMC 
model in a new form that facilitates its computational efficiency. To 
this end, the explicit matrix form of the HFGMC model was 
presented, which allowed an immediate and convenient computer 
implementation of the model. It was shown that a dramatic reduction 
of the computational effort can be achieved by performing a conden- 
sation procedure that, for the doubly periodic HFGMC, reduces the 
dimension of the condensed system of equation to 6 NpN y x 6NpN y , 
just like the reformulation discussed previously. 

III. Progressive Damage Model 

A progressive damage model is sought to function within the 
HFGMC micromechanics model at the level of the fiber/matrix 
constituents. The two significant advances beyond the damage 


models previously employed within HFGMC (cf., Moncada et al. 
[ 6 ]) are 1 ) multiaxiality, such that the response in normal tension, 
normal compression, and shear may be independent, and 2 ) progres- 
sion, such that the damage can occur gradually rather than 
instantaneously. It is also desirable to accomplish these goals as 
simply as possible and while requiring as few model parameters as 
possible. As such, an irreversible elastic damage model is proposed 
that attributes all material nonlinearity to a reduction in material 
stiffness properties. 

The elastic continuum damage model tracks the damage in a 
material element through six scalar damage variables: Df and Df, 
i = 1 , 2, 3. Each damage variable is associated with damage oriented 
normal to the Cartesian coordinate direction indicated by the 
subscript, and separate damage variables track tensile T and 
compressive C damage. While these damage variables may be 
written in vector form, the model is not a vector damage model but 
rather a multiscalar approach, as the damage parameters do not 
necessarily adhere to first-order tensorial mathematics. 

It is assumed that the damage variables modify the engineering 
material properties rather than the stresses or strains, and those in the 
undamaged state, Df = Df = 0, while the completely damaged 
(final failure) state corresponds to Df = 1 or Df = 1. The 
dependence of the engineering material properties on the damage 
variables is taken as 


Ei 


where. 


dil?}-. 

E 2 

= d 2 E° 2 ; E 3 -- 

d 3 E° 3 ; 


V 12 = 

div° l2 - 

V 21 = d 2 v° 2l - 



v 13 = 

d t v 0 l3 - 

v 3 i = d 3 v° 31 ; 



v 23 = 

d 2 v 23 ; 

NJ 

II 

WO 


( 1 ) 


H 

1 -bfDf a u > 

1 - bfDf Ojj < 

0 

0 

( 2 ) 


CT, ; , i= 1, 2, 3 are the normal stress components (no summation), E}, 
£?, and E 3 are the initial material Young’s moduli, Vj 2 , u 21 , v ° l3 , V 31 , 
v 23 , and are the initial material Poisson ratios, and the corre- 
sponding quantities without the 0 superscript are the current 
(damaged) material properties. The bf and bf terms (no summation) 
are constants that control the dependence of the damaged moduli on 
the appropriate damage variables and thus enable scaling of this 
dependence between tension and compression. The default value of 
these b terms is one, and the minimum allowable value of dj is zero, 
but in practice, a small positive minimum value of dj (such as 0 . 000 1 ) 
is employed to promote stability of model simulations. Note that, 
because the Young’s moduli are affected independently by the 
damage variables, a constituent material that is initially isotropic 
becomes orthotropic upon damage initiation. 

It is assumed that the shear moduli dependence on the damage 
variables is given by a linear combination of the appropriate two 
damage variables: 


G 23 = (l-b 42 D^b 43 D' 3 )G° 23 - 
G n = (l-b 5l D’-b 53 Dl)G° l3 ; 

G u = (l-b 61 D\-b 62 DDG° 12 (3) 


where 




Df CT„ > 0 
c ", <0 


( 4 ) 


and GSJj, G° l3 , and G | 2 are the initial material shear moduli. The b 
terms again enable scaling of the damage variable dependence of the 
material properties, and the shear related b terms in Eq. (3) have 
default values of 0.5. The factors that modify the initial shear moduli 
in Eq. (3) (e.g., 1 — b 42 D 2 — b 43 D 3 ), have a minimum allowable 
value of zero, but in practice, a small positive minimum value (such 
as 0 . 0001 ) is employed to promote stability. 
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It is noted that the dependence of the damaged Poisson ratios on 
the damage variables [Eq. ( 1)] preserves the symmetry of the material 
compliance matrix. For example, denoting the compliance matrix 
components by Sy : 


Sn 



v 12 

£7 


(5) 


Substituting for the material properties using Eq. (1_), 


“2^21 

d 2 E° 


“1^12 _ cO _ cO 

' d t E° ~ 12 ~ 21 


(6) 


and, as shown, since the initial compliance matrix is symmetric, it is 
clear that the damaged compliance matrix will remain symmetric as 
well. In addition, S I; = S°j, i j, so the damage only affects the 
diagonal terms of the compliance matrix. Finally, this form of 
damage for the Poisson ratios is mechanistically correct. That is, if 
one considers a typical Poisson mechanical test where, for example, 
a n ^ 0, and all other stress components are equal to zero, we have 
by definition 


where e ;■ are the strain components. Thus, considering the case of 
damage only in the x 1 direction (d { — »• 0), one would expect a small 
g 22 and a large g u , leading to a small v l2 . According to Eqs. (J_) and 
(2), it is clear that the preceding case corresponds to a low value of d t 
and, thus, a small v 12 . Now, considering a constrained mechanical 
test where, for example, a n ^ 0 is applied, s 22 = 0 is constrained, 
and a 22 0 with all other stress components equal to zero, we have 

T S 22 (r 22 — s 22 — 0 (8) 


Thus, substituting using Eqs. (1_) and (5), 


a 22 — — 



E 2 v 12 

E l 




d 2 E°d l v° [ . 


E° v 


■* 2 ^ 2 U 1' , 12 _ j ^2 

■°ii — a 


d i£? 


E° 


(9) 


Therefore, in the case of damage in the x 2 direction (d 2 — »• 0), we see 
that the material loses the ability to accumulate normal stress in the x 2 
direction, regardless of the damage in the x { direction. 


A. Damage Initiation 

The damage initiation criterion is taken as a three-dimensional 
extended version of the strain-based Hashin [ 18 ] criterion. The 
following damage strains are defined. 


are the tensile and compressive damage initiation strains in the 
x 3 -coordinate direction, and Q e ,R e , and S e are the damage initiation 
engineering shear strains associated with y 2 3 , y 13 , and y l2 , 
respectively. Thus, according to Eq. (10), damage will initiate when 
ef , e 2 , or gf exceed one. 

B. Damage Evolution 

The fundamental tenet of the proposed irreversible elastic damage 
model is that the damage evolution is controlled by the damage 
strains, as defined in Eq. (10). That is, while the damage evolution 
controlled by these damage strains can be fit to reproduce uniaxial 
material behavior, it is an assumption of the proposed theory that the 
generalization to multiaxiality can be represented by the three 
damage strains defined in Eq. (10). This is, in many ways, analogous 
to the role of the second stress invariant in von Mises plasticity 
theory, which defines the multiaxial initiation (yield) condition and 
controls the evolution (flow). 

Assuming that the loading imposed on the damaging material is 
done incrementally, as indicated by the dots in Fig. 2, it is possible to 
determine the damage variable increments from the damage strain 
components and damage strain increments given the instantaneous 
slope, kf, of the material’s uniaxial stress-strain response curve. 
Thus, the uniaxial material stress-strain response curve can be used 
to characterize the damage evolution through the damage initiation 
strain and the instantaneous slope. In the uniaxial normal loading 
case with loading in the x f direction, the damage strain reduces to 
gf = g a/X e , and a corresponding damage stress, of\ which is 
linearly related to the damage strain by the damaged modulus, Ef, 
can be identified. Figure 2 graphically shows the relationship for the 
case normal uniaxial loading. 

Because the loading is uniaxial, the modulus at any loading 
increment point is related to the initial modulus E° by the appropriate 
damage variable Dy. 


£? = (! — AOE? (12) 

The incremental change in the modulus can then be determined as 
dEf = -dD t E° (13) 


and the damage variable increment is then given by 


dDj = 


dEf 


(14) 


The damage stress is defined via a linear relationship with the damage 
strain that, for the uniaxial loading, is given by 


where 

X.= 




On > 0 


o-n < O’ 

^ = - 

( z T 

O 33 > 0 

N 

II 

N 

^33 0 


yj a 22 > 0 
Yf o 2 2 < 0 ’ 


(ID 


g u ,i = 1, 2, 3, are the normal stress components (no summation), y 23 , 
y 13 , and y I2 are the engineering shear strains, Xf and Xf are the 
tensile and compressive damage initiation strains in the x { - 
coordinate direction, Yj and Yf are the tensile and compressive 
damage initiation strains in the x 2 -coordinate direction, Zj and Zf 


Efs? 


(15) 



Fig. 2 An incremental elastic damage stress vs damage strain curve. 
Note that X = X e , Y s , or Z,. 
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Note that, even in the uniaxial case, Poisson effects will give rise to 
damage strains in the remaining two normal directions according to 
Eq. (10). However, these will be much smaller than the damage strain 
in the applied loading direction and typically will not exceed one 
(which would indicate damage initiation in these directions) before 
final failure in the applied loading direction. 

From Eq. (15), the incremental change in damage stress is given by 

do? = E?ds? + efdE? (16) 


The instantaneous slope of the damage stress vs damage strain curve, 
which is equivalent to the instantaneous slope of the uniaxial stress- 
strain curve, kj, is given by 


do? 

def 


d FP 

= E? + s?^=k, 


(17) 


Isolating the incremental change in the damaged modulus gives 


dE? = (k, - E?) 



(18) 


and substituting Eq. (18) into Eq. (14) yields 


dDf = 


(E? - k\ ds? 
{ E° )l? 


(19) 


Substituting Eq. (12) into Eq. ( 19 ) and letting k' 
the equation for the damage variable evolution: 


dD l = (1 - D, - £') 



kJE^ provides 


( 20 ) 


Thus, given the fact that the damage initiates when £>, = 0, along 
with the knowledge of the shape (i.e., instantaneous slope) of the 
postinitiation uniaxial stress-strain curve, the incremental damage 
evolution can be determined from Eq. (20). A uniaxial shear stress- 
strain curve could also be used and, if both normal and shear curves 
are available, they can be used to estimate the value of the b terms in 
Eqs. (2) and (3). Otherwise, the default values may be assumed. 


C. Unloading, Reloading, and Load Reversals 

It is assumed that, upon unloading, the material immediately 
ceases to damage. That is, the damage strain state must be on a 
damage surface in order for further damage to occur. As such, the 
maximum previous value of each damage strain must be tracked and 
further damage only permitted when the previous maximum is 
exceeded. Because the model is elastic, unloading will always occur 
linearly to the zero stress-strain state. It is assumed that entirely 
different damage mechanisms are at work in tension and in compres- 
sion. For instance, tensile damage might be due to microcracking, 


while in compression, it might be due to a microbuckling 
mechanism. This implies that both tensile and compressive damage 
strains must be tracked independently, with each causing the 
evolution of its corresponding tensile or compressive damage 
increment, according to Eq. (20). The previous maximum of the 
damage strains in both tension and compression must also be tracked. 
The result of these assumptions is that tensile damage will not affect 
compressive damage, and vice versa. Thus, if the material is 
subjected to a load reversal, the applicable damage variable will 
switch instantaneously, as will the apparent stiffness of the material 
(assuming damage has initiated). For example, if the material is first 
damaged in tension, unloaded, and then subjected to compression, 
the apparent stiffness will change instantaneously to the initial 
stiffness when the material goes into compression. 

The behavior described previously is illustrated in Fig. 3a. The 
material is first loaded in normal tension. The behavior is linear until 
damage initiates at a strain of 0.0125. The postdamage slope is 
negative, and the stress decreases in step 2 as the strain increases. In 
step 3 , the material is unloaded, and the response is linear to the origin 
with a reduced (damaged) stiffness. Upon crossing the origin, the 
material is placed in compression in step 4 and, since the material is 
as yet undamaged in compression, the material returns to its original 
stiffness. Compressive damage initiates at a strain of —0.0287, and 
the compressive postdamage slope in step 5 is positive. The material 
is unloaded in compression in step 6, and the response is again to the 
origin along a path with a reduced (damaged) stiffness. Upon 
reentering the tensile regime in step 7, the material retraces its 
previous tensile unloading path, as the material has already been 
damaged in tension. Further damage does not occur until the previ- 
ously attained maximum strain has been reached, at which point, in 
step 8, the material continues to damage. Clearly, vastly different 
nonlinear material tensile and compressive behavior can be accom- 
modated by the proposed damage model. 

In contrast, example material shear behavior is shown in Fig. 3b. 
Steps 1-3 are similar, but it is clear that in shear (as there is no 
distinction between tension and compression) there is no change in 
slope when traversing from positive to negative stress. Further, 
additional damage in negative shear begins to accumulate when the 
previous maximum positive shear strain (magnitude) is reached. 
When reentering the tensile regime, there is again no change in slope, 
and the previous maximum negative shear strain (magnitude) is 
exceeded before additional damage can occur. Note that Figs. 3a and 
3b are based on the epoxy resin material properties for which the 
characterization is discussed in Sec. IV. 


D. Final Failure 

Final failure is considered to be the state at which the material can 
no longer support one or more components of stress in tension, 




a) b) 


Fig. 3 


Monolithic material described by the proposed damage model: a) cyclic tensile and compressive behavior and b) cyclic shear behavior. 
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compression, or shear to any significant degree. At this point, the 
effective stiffness of the material, rather than being determined 
according to Eqs. (J_) or (3), will be set to a very low value (i.e., 
0.0001 times its initial value). One way that the material can reach 
final failure is by one of its damage variables evolving to the point that 
one of the stiffnesses reaches zero. It is desirable, however, to enable 
the material to reach final failure suddenly before a stiffness is 
evolving all the way to zero. 

Separate final-failure criteria are proposed for tension and 
compression. In tension, a strain energy release rate criterion is used. 
For a body with a macroscopic crack, the strain energy release rate G 
is the amount of energy released in advancing the crack by a distance 
da normalized by the new crack area. It is typically written as 


where t is the depth of the cracked body (thus, tda is the new crack 
area). In adapting the concept of strain energy release rate to 
continuum damage mechanics, we work with the strain energy 
density released W rather than the total energy released W s Thus, we 
use the relation 

dW s = VdW (22) 

where V is the volume of the material that is being damaged, and the 
increment in released strain energy density is shown graphically in 
Fig. 4. By assuming a piecewise linear response, as in Fig. 4, the 
released strain energy density can be calculated as 


Table 1 Association between loading modes and stress 
components for damage associated with the 
coordinate direction, x, 


Coordinate direction 


Stress component associated 
with loading mode 

i (damage orientation) 

j 

k 

i ii m 

i 

2 

3 

^11 r 12 T 13 

2 

3 

1 

°22 T 23 T 12 

3 

1 

2 

°33 T 13 T 23 


summation), mode II damage is caused by r * •, and mode III damage is 
caused by x ik (see Table 1_). The corresponding mode-specific strain 
energy density increments are then given by the increments in 
released energy density, as shown in Fig. 5, associated with the 
(J u — s u , Ty - Yjj, and x ik — y ik response curves. 

Because, as indicated in Table T each shear stress component is 
associated with a loading mode for damage in two directions (i.e., x l2 
is associated with mode II loading for damage in the x, direction as 
well as mode III loading for damage in the x 2 direction), it is 
necessary to partition the shear strain energy density increments 
appropriately. Toward this end, denoting the stress and strain in 
vector form, 

a = [a n a 22 a 33 r 23 t 13 t 12 ] 7 ; 

e = [fill ®22 e 33 f23 Xl3 Yllf (24) 


dW = |{ff(e + ds) — e(a + da)} (23) 

It is possible to separate the strain energy release rates associated 
with each tensile damage mode, where again an analogy with a 
macroscopically cracked body is used. The three modes of loading 
for such a macroscopic crack are denoted as I, II, and III, and they are 
associated with opening, in-plane shear, and out-of-plane shear, 
respectively [26], Figure 5 illustrates the three modes of loading as 
they apply to tensile damage in the x, direction, which would be 
associated with the damage variable Df, where i, j, k = 1, 2, 3, and 
i ^ j ^ k. In this scenario, mode I damage is caused by a u (no 



Fig. 4 Graphical representation of the released strain energy density 
increment. 



Mode I Mode II 


Mode III 


Fig. 5 Three loading modes that affect damage in the x, direction. 


and denoting the components of the stress and strain vectors as a„ and 
e„, respectively, with n = 1, .... 6, the strain energy density release 
rate associated with each response curve can be denoted as dW n . The 
following partitions can be made according to damage direction in 
order to arrive at mode-specific strain energy density release rates: 


dW\ = dWy, 


dW],= 


dD , 


dD\ + dD 2 


-dWf. 


dW)„ = 


dD, 


dWji = 


dD y + dD-y 
dD 2 


dW , 


dD 2 — dD 


■dW, 


dWn, = 


dWj = dW 2 ; 
dD 2 


dDy + dD 2 


dWf.: 


dW, = dW 3 ; dWn = 


dD 3 

dDy + dD 3 


dW<: 


dwi„ = 


dD , 


dD~, + dD 


-dWi 


(25) 


where dW' M is the mode-specific increment in released strain energy 
density for damage in the x ; direction and M is the mode (I, II, or III). 

The mode-specific strain energy release rates associated with 
damage in the x ; direction can be found using Eqs. (21), (22), and 
(25) and are given by 


C'u = 


dW' M 

hdaV 


-V = - 


dW‘ M 

hda\ 


111 dW ‘ M ll 
ijk dd< ij 


(26) 


where, as shown in Fig. 5, the lengths of the material in the three 
directions are denoted by /„ lj, and l k , and da' M is a mode-specific 
increment in effective crack length, which can be calculated by 
scaling the length of the material in the damage direction by the 
increment in stiffness reduction due to the damage variable 
increment. That is, 

d a j ~ l 2 b y yd l) ] , da H — l 2 b^ydDy, dci k22 — l*>b 3 ydD j, 

dd] = l 3 b 22 dD 2 ; dcr u = l 3 b 62 dD 2 , dcr UI = l 3 b K dD 2 ; 

da] = lyb 33 dD 3 , da 3 n = l t b A 3 dD 3 \ da] u = Iyb 53 dD 3 (27) 


Substituting Eq. ( 27 ) into Eq. (26) yields 
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h 

dW] 


h dW), 


h 


b T u 

d£>, ’ 

u ii 

b&i dD\ 

U III 

bsi 

dO, 

h 

dW] 

rz 2 

l 2 d W], _ 

rz 2 

h 

dW] n 

b 2 2 

dD 2 ’ 

u ii 

b 42 dZ) 2 

u iii 

b(, 2 

dD 2 

h 

d W) 

rz 3 

/ 3 dW 3 ,, 


_ h 

d W)„ 

b\ 3 

d D 3 ’ 

u ii 

Z?53 dD 3 

U III 

b<n 

d£>, 


which enables the calculation of the instantaneous mode-specific 
strain energy release rate for each damage direction (indicated by the 
superscript). 

The strain energy release rates from Eq. ( 28 ) can then be used in 
final-failure criteria by comparing them to the material critical strain 
energy release rates. For example, one can use the maximum-strain 
energy release rate criterionm in which final failure occurs when any 
of the following are satisfied: 

G' m = G c m (29) 


where is the mode-specific critical strain energy release rate. 
Alternatively, the power-law criterion is given by 



where a is a material parameter. Finally, the Benzeggagh-Kenane 
(B-K) [27] criterion can be employed, which is applicable only when 
Gfj = G^n- The B-K criterion can be expressed as 

G‘ + G' n + G' ln = G c (31) 


with 


G c = G i + (Gj, — Gf ) 


( G), + G)„ y 

\G‘, + G\, + G‘J 


(32) 


where is a material parameter. 

In compression, the failure mechanism is assumed not to be related 
to cracking, and thus a strain energy release rate based final-failure 
criterion is not applicable. Instead, a criterion based on the total 
dissipated energy associated with a damage direction is employed. 
This criterion can be expressed as 

(W\ + W i II + W i m )V=WC (33) 

where W!r is a critical strain energy and W' M is the mode-specific 
strain energy release rate, which can be determined by integrating the 
mode-specific strain energy release rate increments. 


IV. Results and Discussion 

To demonstrate the capabilities of the HFGMC micromechanics 
model, including the proposed progressive damage model, a 
commonly modeled composite system has been chosen. This is the 
£-glass/MY750/HY917/DY063 epoxy resin composite that was 
included in the WWFE [JJ. As is often the case, while nominal 
constituent (fiber/matrix) elastic and strength properties were 
provided in the WWFE. nonlinear stress-strain curves of the 
constituents, as are desirable for nonlinear micromechanics model 
input, were not provided. As such, the constituent properties were 
obtained based on correlation with composite ply-level nonlinear 
data provided in the WWFE for the £-glass/epoxy composite (i.e., 
the constituent properties were backed out). 

All of the results presented consider a 7 x 7 subcell RUC with a 
fiber volume fraction of 0.6, as shown in Fig. 6. While HFGMC is 
mesh-dependent, this HFGMC unit cell discretization was shown to 
yield a converged global nonlinear deformation response for a SiC/Ti 
composite by Bednarcyk et al. [ 28 ] Note that the doubly periodic 
RUC shown is infinitely long in the Aj direction and, thus, represents 
a special case of the triply periodic version of HFGMC described 
previously. It should also be noted that, for nonlinear analysis, 
HFGMC tracks the local fields within each subcell at a number of 


























Fiber 


























- Matrix 


Fig. 6 7x7 RUC used to represent the composite material. 


integration points. Three integration points were used for this 
purpose in both the x 2 and x 3 directions, for a total of nine integration 
points per subcell. Obviously, damage will initiate and evolve 
differently at each of these nine points per subcell; thus, an average of 
the damage at these nine points was used to determine the subcell 
damaged elastic properties. 

The simulations reported next were executed on a standard dual- 
core Windows workstation. The execution times for generating the 
stress-strain response of the neat epoxy was a fraction of a second 
(~0. 1 s), whereas, the stress-strain curves for the composite (using 
the 7 x 7 unit cell) required approximately 3 s of execution time. 
Generation of the composite failure envelopes was accomplished by 
repeatedly simulating the stress-strain response along radial paths 
with a given angular increment in the given biaxial stress space under 
stress control. That is, to generate a simulated failure envelope in 
a, | — a 22 stress space, a simulation was first run by monotonically 
increasing a n with all other stress components equal to zero until 
final failure was predicted. Then, a biaxial stress-strain simulation 
(monotonically increasing both a x , and a 22 simultaneously) was run, 
where tan _1 (a 2 2/°'ii) = 5° (for example) until final failure. This 
process continued until the entire envelope was generated through 
the complete 360° in the stress space. Because the failure envelope 
generation is automated within the HFGMC code, these composite 
failure envelope simulations typically require only several minutes 
(~4 min) of execution time. 

Figures 7a-7c show the correlation of the micromechanics model 
with the composite ply-level experimental data, provided by Hinton 
et al. [JJ, used to determine the constituent properties. The 
progressive damage model assumes that all material nonlinearity is 
due to stiffness reduction rather than inelastic deformation. 
Unloading behavior of the composite, were it available, could be 
examined to assess the validity of this assumption. The determined 
constituent data are summarized in Tables 2 and 3. While the resin 
matrix is subject to the proposed progressive damage model, the fiber 
is treated as isotropic and linear elastic until failure. A maximum 
stress criterion is used for the fiber, based only on tensile and 
compressive longitudinal stresses. Once the fiber fails, all of its 
elastic stiffness properties are instantaneously reduced to a very low 
value (i.e., D' are set to 0.9999, see Eqs. (1-4)). It is noted that the 
fiber strengths determined via correlation in Table 2 are slightly lower 
than the nominal values of 2 1 50 and 1450 MPa given by Hinton et al. 
[JJ Further, in the model longitudinal composite tensile response, 
shown in Fig. 7a, it is noted that some nonlinearity due to matrix 
damage is present before failure. 

As stated, the matrix is assumed to be subject to the proposed 
progressive damage model. An exponential form was assumed for 
the instantaneous modified slope [see Eq. (20)] of the matrix material 
stress-strain curve after damage initiation: 

k' = Ae~ s ? /B (34) 

where the constants A and B are permitted to take on different values 
in tension T and compression C. As indicated in Table 3, these 
constants were determined via correlation with the ply-level 
nonlinear composite response curves provided by Hinton et al. [JJ 
For final failure of a subcell in tension, the maximum-strain energy 



a n (MPa) 
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a) b) 



e 12 


c) 

Fig. 7 Model correlation with experiment for a) the longitudinal tensile and compressive response, b) transverse tensile and compressive response, and 
c) the longitudinal shear response of /i-glass/epoxy. 


Table 2 Properties of Silenka /J-glass fiber 


Property 

Symbol 

Units 

Source 

Value 

Elastic modulus 

E 

GPa 

Hinton et al. [1] tabular value 

74 

Poisson’s ratio 

V 

— 

Hinton et al. [1] tabular value 

0.2 

Longitudinal tensile strength 

x T 

MPa 

Correlation (Fig. 7a) 

2110 

Longitudinal compressive strength 

x c 

MPa 

Correlation (Fig. 7a) 

1290 


Table 3 Properties of MY750/HY917/DY063 epoxy 


Property 

Symbol 

Units 

Source 

Value 

Elastic modulus 

E 

GPa 

Correlation (Fig. 7b) 

3.7 

Poisson’s ratio 

V 

— 

Hinton et al. [1] tabular value 

0.35 

Tensile damage initiation strain 

Xj 

— 

Correlation (Fig. 7b) 

0.0125 

Compressive damage initiation strain 

X e c 

— 

Correlation (Fig. 7b) 

0.0287 

Engineering shear damage initiation strain 

S e 

— 

Correlation (Fig. 7c) 

0.0443 

Postdamage slope parameter, Eq. (34) 

a t 

— 

Correlation (Fig. 7c) 

0.7 

Postdamage slope parameter, Eq. (34) 

A c 

— 

Correlation (Fig. 7b) 

2.0 

Postdamage slope parameter, Eq. (34) 

b t 

— 

Correlation (Fig. 7c) 

0.82 

Postdamage slope parameter, Eq. (34) 

B c 

— 

Correlation (Fig. 7b) 

0.96 

Tensile scaling parameter, Eq. (2) 

bl 

— 

Correlation (Figs. 7b and 7c) 

1.32 

Compressive scaling parameter, Eq. (2) 

bf, 

— 

Correlation (Figs. 7b and 7c) 

1.32 

Shear scaling parameters, Eq. (3) 

bn, b 5i , b 6i 

— 

Default value 

0.5 

Mode I critical strain energy release rate 

G? 

J/m 2 

Correlation (Fig. 7b) 

800 

Modes II and III critical strain energy release rates 

Gf,(=Gf„) 

J/m 2 

Correlation (Fig. 7c) 

2400 

Critical compressive strain energy 

wf 

j 

Correlation (Fig. 7b) 

1.86 x 10~ 6 

Material length 

l, 

m 

Correlation (Figs. 7b and 7c) 

9 x 10- 5 
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release rate criterion [Eq. (29)] was used. It should also be noted, in 
Table 3, that the elastic modulus value employed for the matrix is 
higher than the 3.35 GPa nominal value provided by Hinton et al. [1] 
Further, the tensile and compressive scaling parameters 
were chosen to scale the effect of normal damage with respect to 
shear damage. The default shear scaling parameters (b v , b 5i , b 6i ) 
were employed. These parameters provide the damage model with 
the capability to behave differently in shear vs normal deformation, 
despite the fact that identical A and B parameters control the damage 
evolution in both cases. Hence, one can correlate the normal response 
using default scaling parameters and then choose the shear scaling 
parameters to enable capture of the shear response, or vice versa (as 
done here). As the damage model is presented herein, the scaling 
parameters can also be assigned different values per stiffness 
component, but this would result in anisotropic progressive damage 
behavior of an initially undamaged material and is thus not desirable 
in the present application to the neat resin matrix. Finally [from 
Eqs. (28) and (33)] it is clear that the critical strain energy release 
rates and the critical compressive strain energy are scaled by the 
material lengths. Therefore, these parameters are not independent 
and, by choosing a different material length, different values of the 
critical strain energy release rates and critical compressive strain 
energy would result from the correlation with the ply-level composite 
data. 

Figure 8 shows the monolithic resin matrix behavior that results 
from the parameters given in Table 3. Most noteworthy is the 
extremely brittle behavior in tension compared with the ductile- 
appearing behavior in compression and shear. This brittle behavior is 
a direct result of the brittle transverse tensile composite behavior 
observed in Fig. 7b that was used to characterize the matrix tensile 
damage/failure response. Obviously, it would not be possible to 
capture these vast differences in the character of the damage/failure 
behavior among tension, compression, and shear via a simple failure 
criterion that treats the epoxy as linear to failure, as shown in [6], 
Similarly, a nonlinear model similar to plasticity, which characterizes 
nonlinearity through a single scalar measure (i.e., von Mises stress), 
would miss the brittle character in tension if it were to capture the 
more ductile-appearing character in compression and shear. Thus, in 
the context of micromechanics, the need for a damage model that can 
handle material nonlinearity and also separate the damage modes is 
clear, as is offered by the present model. 

Figures 9-1 1 provide the model damage initiation and final-failure 
envelopes for the monolithic (neat) epoxy resin, given the material 
properties in Table 3. The distance between the two envelopes in each 
figure indicates the degree to which the predicted failure is 
progressive vs sudden. For instance, in Fig. 9, it is clear that, in the 
tension-tension regime (quadrant I), the failure is very brittle, with 
virtually no progression of damage between initiation and final 
failure. In the compression-compression regime (quadrant III), 
however, there is a great deal of progression before final failure. 



Fig. 8 Monolithic (neat) epoxy resin behavior simulated via the 
proposed damage model, given the material properties shown in Table 3. 



Fig. 9 Model damage initiation and final-failure envelopes for the 
monolithic (neat) epoxy in the biaxial normal stress plane. 


Figure H) provides the normal-shear failure envelopes for the case 
in which the normal and shear components are interactive. That is, 
examining Eq. (10), it is clear that the £ n and s 12 strain components 
are interactive, in that they both contribute to the damage strain ef. 
This interaction between the 11 -normal and 12-shear behavior is 
manifested as the elliptical character of the surfaces in Fig. 10. 
Figure J_1 , on the other hand, provides the normal-shear envelopes for 
the case in which the components do not interact. Again examining 
Eq. (10), it is clear that , contributes only to damage strain ef , while 
g 23 contributes to ef and ef ; thus, these components do not interact. 
The result is the more rectangular character of the surfaces in Fig. 11 . 
Thus, the material, which was initially isotropic, is highly anisotropic 
in terms of its multiaxial damage progression and failure response. It 
is noted that the damage surfaces in Fig. U_ do display some 
interaction effects because the surfaces are generated in stress space, 
so application of ay t with a 22 = 0 results in nonzero e 22 and g 33 due to 
Poisson effects. 

Another important feature of the implementation of the proposed 
progressive damage model is illustrated in Fig. 10. When the normal 
stress component transitions from tensile to compressive (quadrant I 
to quadrant II transition), there is no jump in the final-failure 
envelope as one might expect (as the compressive damage 



Fig. 10 Model damage initiation and final-failure envelopes for 
the monolithic (neat) epoxy in the normal-shear stress plane with 
interaction. 
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Damage Initiation 

— Final Failure 



Fig. 11 Model damage initiation and final-failure envelopes for the 
monolithic (neat) epoxy in the normal-shear stress plane without 
interaction. 



Final Failure 

— Damage Initiation 
□ Experiment (Characterization) 


o,, (MPa) 


Fig. 13 Model damage initiation and final-failure envelopes for 
unidirectional E-glass/epoxy in the longitudinal tension: in-plane shear 
stress plane. 


parameters in Table 3 become active). Rather, the transition is 
gradual. This feature was incorporated by basing the transition in 
properties on the degree to which the damage strain component is 
shear dominated. That is, when a damage strain is completely shear 
dominated (i.e., on the y axis of Fig. _10), the tensile parameters are 
used. If the stress state is compressive, the parameters are inter- 
polated between the tensile and compressive values until the damage 
strain is no longer dominated by shear, at which point they become 
the compressive values. 

Figures 12-14 show model composite damage initiation and final- 
failure envelopes. The experimental data are taken from Hinton et al. 
[ 1 ] The experimental data on the axes are the data used to characterize 
the damage model for the matrix material. The lack of agreement in 
quadrant IV indicates that a different final-failure criterion, such as 
the power-law or B-K criteria, should be employed. The interaction 
of strain energy release rates provided by these criteria may improve 
the agreement. In Fig. 14, the noncharacterization experimental data 
are actually for a 0.60 volume fraction £-glass/LY556/HT907/ 
DY063 epoxy composite [1], for which the properties are similar to 
the £-glass/MY750/HY917/DY063 epoxy composite, for which the 
materials employed in the simulation were characterized. As shown, 
the agreement is quite good, most notably in the ability of the model 
to capture the observed strengthening in quadrant II. 

As a final example, the progressive damage model has been 
characterized based on the in situ matrix shear response measured 
and presented by Ng et al. [29] The specific type of fiber and 
composition of the polymer matrix were not disclosed by these 
authors for proprietary reasons. Using a novel technique, Ng et al. 
[29] have detennined the average nonlinear in-plane shear behavior 



Final Failure 

— Damage Initiation 
□ Experiment (Characterization) 
o Experiment 


a lt (MPa) 


Fig. 12 Model damage initiation and final-failure envelopes for 
unidirectional /'.'-glass/epoxy in the biaxial normal in-plane stress plane. 



Fig. 14 Model damage initiation and final-failure envelopes for 
unidirectional £-glass/epoxy in the transverse tension: in-plane shear 
stress plane. 


90 n 



Y.2 

Fig. IS Model correlation with polymer matrix shear response 
measured by Ng et al. [29]. 
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Table 4 Properties of the composite constituents considered by Ng et al. [29] 


Property 

Symbol 

Units 

Source 

Value 

Fiber axial elastic modulus 

E n 

GPa 

Ng et al. [29] tabular value 

310 

Fiber transverse elastic modulus 

£*22 

GPa 

Ng et al. [29] tabular value 

20.5 

Fiber major Poisson’s Ratio 

v 12 

— 

Ng et al. [291 tabular value 

0.28 

Fiber minor Poisson’s Ratio 

y 23 

— 

Ng et al. [29] tabular value 

0.35 

Fiber axial shear modulus 

G \2 

GPa 

Ng et al. [29] tabular value 

32.8 

Matrix elastic modulus 

E 

GPa 

Ng et al. [29] 

3.078 

Matrix Poisson’s ratio 

V 

— 

Ng et al. T291 

0.35 

Matrix engineering shear damage initiation strain 

s e 

— 

Correlation (Fig. 15) 

0.02 

Matrix postdamage slope parameter, Eq. (34) 

A T 

— 

Correlation (Fig. 15) 

1.1 

Matrix postdamage slope parameter, Eq. (34) 

B r 

— 

Correlation (Fig. 15) 

3.3 

Matrix shear scaling parameters, Eq. (3) 

b«h b 5i , b 6i 

— 

Default value 

0.5 

Matrix modes II and III critical strain energy 
release rates 

GF,(=Gf,/) 

J/m 2 

Correlation (Fig. 15) 

1240 

Matrix material length 

h 

m 

Correlation (Fig. 15) 

9 x 10~ 5 


of the composite at the ply level, as well as that of the matrix. Using 
the four matrix stress-strain curves given by these authors, the 
correlation of the progressive damage model with the measured data 
is presented in Fig. 15. The corresponding material parameters for the 
polymer matrix are given in Table 4, along with the fiber material 
properties provided by Ng et al. [29]. It is evident from Fig. _F5 that 
the simple exponential form of the instantaneous slope given in 
Eq. ( 34 ) captures the general character of the matrix nonlinear shear 
response reasonably well, but it is unable to precisely capture the 
early stiffness decrease of the matrix as well as the saturation 
behavior close to failure. A different functional form for the 
instantaneous stiffness could provide more precise correlation in this 
regard. Note that the failure point for the model was chosen as the 
average ultimate shear strength of the four measured shear curves 
(80.4 MPa). This point corresponded to a mode II/III strain energy 
release rate of 1240 J /m 2 . Note also that it was assumed for both the 
matrix and the composite that the only contribution to failure was 
matrix shear damage and failure; thus, some model data that were 
given in the previous examples are not needed in the present 
simulations. 

A comparison of the predicted composite in-plane shear response 
with the data provided by Ng et al. [29] is given in Fig. Uj. As before, 
the HFGMC RUC shown in Fig. 6 was used, with nine integration 
points per subcell. As was the case for the correlated matrix shear 
response (Fig. 1_5), the model does a reasonably good job of 
predicting the composite response. The slight overprediction of the 
measured curves near damage initiation and the inability to capture 
well the saturation behavior is also similar to what was observed in 
the model characterization for the matrix. The failure prediction for 
the composite was taken as the point at which the average matrix 
strain energy release rate reached 1240 J/m 2 . This approach was 


90 



Fig. 16 Comparison of the predicted Vf = 68 % polymer matrix 
composite in-plane shear response with that measured by Ng et al. [291. 


used because the critical strain energy release rate was determined 
from the matrix shear response curves given by Ng et al. [29], which 
provide the average in situ matrix behavior. If the failure prediction 
was based on failure of the first subcell, it would have occurred at a 
very low stress-and-strain value (r 12 = 65.3 MPa, y l2 = 0.021), 
along the same model curve shown in Fig. _16. The average matrix 
strain energy release rate approach overpredicted both the ultimate 
shear strength and shear strain to failure of the composite. However, 
this may be due to the fact that the matrix and composite saturation 
behavior have not been well captured by the exponential form for the 
instantaneous slope used in the simulations. Aspects of the 
techniques employed by Ng et al. [29] to determine the matrix shear 
response could contribute to some of the discrepancy as well. 

V. Conclusions 

The HFGMC micromechanics model has been extended to 
include progressive damage at the constituent (subcell) level. The 
proposed coupled micromechanics-damage model seeks to include 
multiaxiality, such that the character of the constituent behavior can 
be vastly different in tension, compression, and shear in as simple a 
way as possible. It appears that enabling this distinct behavior based 
on type of loading is beneficial to predicting the damage and failure 
behavior of polymer matrix composites. The proposed damage 
model attributes all material nonlinearity to stiffness reduction and is 
thus less applicable to materials that exhibit plastic deformation. The 
model’s true applicability to a given material can be assessed by 
examining the material’s unloading response to determine the 
proportion of the nonlinear deformation that is permanent. The 
damage model introduces six damage parameters (three for tension 
and three for compression) and associated damage strains based on 
oriented continuum damage within the constituent materials. 
Assuming that the damage strains adequately capture the multiaxial 
nature of the damage, evolution equations for the damage parameters 
can be established based on the postdamage initiation shape of the 
constituent stress-strain response. Final failure of a constituent 
subcell has also been addressed by calculating mode-specific strain 
energy release rates in tension and a total dissipated strain energy in 
compression. The coupled HFGMC-damage model has been 
characterized for an £-glass/epoxy composite based on data in the 
literature that exhibit brittle behavior in transverse tension but 
ductile-appearing (i.e., nonlinear) behavior in transverse compres- 
sion and longitudinal shear. It was shown that the HFGMC-damage 
model can capture this behavior, and neat resin and composite-level 
stress-strain curves, damage initiation envelopes, and final-failure 
envelopes were presented for this material system. A final example 
was given in which the HFGMC-damage model was applied to 
predict the shear response of a composite for which the matrix shear 
response was available in the literature. Future investigation of the 
presented model will focus on extension to multiscale analysis of 
resin matrix composite laminates and structures, the latter of which 
will necessitate linkage of the couple micromechanics-damage 
model with a finite-element solver. 
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